CB-mouthColor

library(RColorBrewer); #brewer.pal(n=8, name = "Dark2") #;display.brewer.pal(n = 8, name = "Dark2")
library(tidyverse); theme_set(theme_classic()) #ggplot2, dplyr, tibble, etc.
library(plotly) #use ggplotly() for interactive plots with scoll-over IDs
library(knitr) #use kable() to make formatted tables
library(kableExtra)
library(Rtsne)
library(glmmTMB)
library(lme4)
library(sjPlot)
library(ggpubr) #ggarrange() figure wraps
library(plyr)
library(ggeffects) #ggpredict
library(rstanarm) #stan models
library(shinystan) #stan model evaluation >launch_shinystan_demo()
library(loo) #loo() to compare fits between bayesian models
library(MuMIn) #dredge, model averaging
library(AICcmodavg) #aictab() for AICc model comparisons
library(coin) #permutation test for ratios with independence_test()
library(janitor) #get_dupes() finds and prints duplicates

Read in mouth color data (mc)

allAges.mc=tibble(read.csv("CB-mouthColor.csv", h=T)) #N=215

#subset with ages 23-27
mc <- allAges.mc %>% filter(between(age,23,27)) #N=133

Create subsets of mc for females (N=78) and males (N=55) only

#females only
female.mc <- mc %>% 
  filter(sex == "F") #N=78
#males only
male.mc <- mc %>% 
  filter(sex == "M") #N=55

Read in nestling morphometrics data (nm)

Nestling morphometrics are being used here to boost sample size for calculating mass/tarsus residuals
Replace zeros with NAs, sort and filter by age (23-27)

nm=tibble(read.csv("nestlingMorph.csv", h=T))
#replace zeros with NAs in morphometrics 
nm[nm==0] <- NA
#remove nestlingMeasurementNumber
nm <- nm[,-1]
#sort by year > nestName > id
nm <- nm %>% arrange(year, nestName ,id)
#filter to include only ages 23-27
nm <- nm %>% filter(between(age, 23, 27)) 

Create subset of nm including only sexed birds (N=503; females=268; males=235)

#subset nm to include only sexed birds and remove NAs from mass & tarsus
sexed.nm <- nm %>% drop_na(sex, tarsus, mass) #N=503
#make sex a factor rather than character
sexed.nm$sex <- as.factor(sexed.nm$sex)

Plot–mass distribution by sex

Only includes sexed individuals (N=503)

sexed.nm %>% 
  ggplot(aes(mass, fill=sex))+
  geom_histogram(binwidth = 10)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))+
  annotate(geom = "text", x=375, y=29, label="N=268")+
  annotate(geom = "text", x=390, y=5, label="N=235")

Plot–tarsus distribution by sex

Only includes sexed individuals (N=503)

sexed.nm %>% 
  ggplot(aes(tarsus, fill=sex))+
  geom_histogram(binwidth = 0.5)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))+
  annotate(geom = "text", x=59, y=28, label="N=268")+
  annotate(geom = "text", x=60, y=5, label="N=235")

Plot–mass by tarsus with sex filter (N=503)

massTarsus.plot <- sexed.nm %>% 
  ggplot(aes(x=tarsus,y=mass,color=sex,label=id))+
  geom_point(alpha=0.5)+
  geom_smooth()+
  scale_color_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))
ggplotly(massTarsus.plot)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Remove outliers (SWB-W BURP03, EN WKAY05, GL JSUP05, YASB WCEM01), now N=499

sexed.nm <- sexed.nm[-c(211, 141, 332, 353),] #N=499

massTarsus.plot <- sexed.nm %>% 
  ggplot(aes(x=tarsus,y=mass,color=sex,label=id))+
  geom_point(alpha=0.5)+
  geom_smooth()+
  scale_color_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))
ggplotly(massTarsus.plot)

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

#females only 
female.nm <- sexed.nm %>% 
  filter(sex == "F") #N=265
#males only
male.nm <- sexed.nm %>% 
  filter(sex == "M") #N=234

Create the female mass~tarsus residual index

femaleIndex.lm <- lm(mass~tarsus, female.nm)
summary(femaleIndex.lm)

## 
## Call:
## lm(formula = mass ~ tarsus, data = female.nm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -107.065  -17.796    2.143   18.200   95.956 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -373.9956    44.1717  -8.467  1.8e-15 ***
## tarsus        12.6523     0.7617  16.610  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 30.89 on 263 degrees of freedom
## Multiple R-squared:  0.512,  Adjusted R-squared:  0.5101 
## F-statistic: 275.9 on 1 and 263 DF,  p-value: < 2.2e-16

Check residual distributions

*Q-Q plot looks a little non-normal, probably seeing the effects of reduced sample size by separating sexes.

plot(femaleIndex.lm)

### Create male mass~tarsus residual index

maleIndex.lm <- lm(mass~tarsus, male.nm)
summary(maleIndex.lm)

## 
## Call:
## lm(formula = mass ~ tarsus, data = male.nm)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -114.779  -24.276    2.291   26.091  114.680 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -334.1798    52.9473  -6.312 1.39e-09 ***
## tarsus        12.0998     0.8885  13.618  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 38.43 on 232 degrees of freedom
## Multiple R-squared:  0.4442, Adjusted R-squared:  0.4419 
## F-statistic: 185.5 on 1 and 232 DF,  p-value: < 2.2e-16

Plot residuals

plot(maleIndex.lm)

Create dataframes containing id and residuals

#male residuals pulled from lm() output
male.nm$resids <- maleIndex.lm$residuals
#female residuals pulled from lm() output
female.nm$resids <- femaleIndex.lm$residuals

#dataframes containing male and female ids and residuals for the join with mc
joinMaleResidual <- data.frame(id = male.nm$id,
                               resids = male.nm$resids)
joinFemaleResidual <- data.frame(id = female.nm$id,
                                 resids = female.nm$resids)

Join residuals with male and female mouth color data by common id

#left join with female mc
female.mc <- left_join(female.mc,joinFemaleResidual, by="id")

#left join with male mc
male.mc <- left_join(male.mc,joinMaleResidual, by="id")

#left join with both sexes (full mc) 
bothSexResid <- data.frame(rbind(joinMaleResidual,joinFemaleResidual))
mc <- left_join(mc,bothSexResid, by="id")

Distribution of residuals for both sexes combined

mc %>%
  ggplot(aes(x=resids, fill=sex))+
  geom_histogram(binwidth = 10)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))

Mass and mass/tarsus indices for modeling

mc$ratio <- mc$mass/mc$tarsus
mc$cubeMass <- mc$mass^(1/3)
mc$cubeRatio <- ((mc$mass^(1/3))/mc$tarsus)*100

Permutation test to assess independence of male and female ratios

reference: https://stats.stackexchange.com/questions/23152/test-for-significant-difference-in-ratios-of-normally-distributed-random-variabl

independence_test(ratio~as.factor(sex), data = mc)

## 
##  Asymptotic General Independence Test
## 
## data:  ratio by as.factor(sex) (F, M)
## Z = -2.6339, p-value = 0.008441
## alternative hypothesis: two.sided

Distribution of mass/tarsus ratio for both sexes combined

mc %>%
  ggplot(aes(x=ratio, fill=sex))+
  geom_histogram(binwidth = 0.4)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))

mc %>%
  ggplot(aes(x=mass, fill=sex))+
  geom_histogram(binwidth = 20)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))

mc %>%
  ggplot(aes(x=cubeMass, fill=sex))+
  geom_histogram(binwidth = 0.1)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))

mc %>%
  ggplot(aes(x=cubeRatio, fill=sex))+
  geom_histogram(binwidth = 0.3)+
  scale_fill_manual(values = c("#1B9E77", "#D73027"),
                    name = "",
                    labels = c("Female","Male"))

SATURATION 1 MODELS

Plot saturation by age

mc %>% 
  ggplot(aes(x=sat1,y=age)) +
  geom_point()+
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Plot saturation by mass

satMass.plot <- mc %>% 
  ggplot(aes(x=sat1,y=mass, label = id)) +
  geom_point()+
  geom_smooth()
ggplotly(satMass.plot) 

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Plot age by mass

ageMass.plot <- mc %>% 
  ggplot(aes(x=age, y=mass, label = id)) +
  geom_point()+
  geom_smooth()+
  geom_boxplot()
ggplotly(ageMass.plot) 

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

## Warning: Continuous x aesthetic -- did you forget aes(group=...)?

Sat1 ~ mass/tarsus residual model

sat1Residual.mdl <- lmer(sat1 ~ age + resids + sex + 
                           bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = mc, REML = FALSE)
summary(sat1Residual.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ age + resids + sex + bodTemp1 + ambTemp1 + timeOut1 +  
##     (1 | family)
##    Data: mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    805.4    828.9   -393.7    787.4       91 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7578 -0.5450  0.1383  0.5430  2.1740 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 108.51   10.417  
##  Residual              94.31    9.711  
## Number of obs: 100, groups:  family, 36
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 126.33319   73.96937   1.708
## age          -3.81740    1.54508  -2.471
## resids        0.08675    0.04658   1.862
## sexM         -2.14906    2.41129  -0.891
## bodTemp1      1.64250    1.86653   0.880
## ambTemp1      1.30666    0.59721   2.188
## timeOut1     -0.01112    0.14501  -0.077
## 
## Correlation of Fixed Effects:
##          (Intr) age    resids sexM   bdTmp1 ambTm1
## age      -0.292                                   
## resids    0.133 -0.242                            
## sexM     -0.140 -0.069  0.106                     
## bodTemp1 -0.836 -0.260  0.020  0.154              
## ambTemp1  0.093  0.202 -0.160  0.039 -0.345       
## timeOut1 -0.027 -0.088  0.043  0.001  0.066 -0.199

Plot residual model

plot_model(sat1Residual.mdl,
           title = "Saturation 1 residual model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Sat1 ~ mass/tarsus ratio model

sat1Ratio.mdl <- lmer(sat1 ~ age + ratio + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = mc, REML = FALSE)
summary(sat1Ratio.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ age + ratio + sex + bodTemp1 + ambTemp1 + timeOut1 + (1 |  
##     family)
##    Data: mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    805.3    828.8   -393.7    787.3       91 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6577 -0.5549  0.1515  0.5213  2.0561 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 99.29    9.964   
##  Residual             97.09    9.853   
## Number of obs: 100, groups:  family, 36
## 
## Fixed effects:
##               Estimate Std. Error t value
## (Intercept) 93.6305265 73.4209470   1.275
## age         -3.6954789  1.4948678  -2.472
## ratio        5.0348153  2.5351722   1.986
## sexM        -3.8538809  2.4927588  -1.546
## bodTemp1     1.5537491  1.8767843   0.828
## ambTemp1     1.4185318  0.5787127   2.451
## timeOut1     0.0005294  0.1435943   0.004
## 
## Correlation of Fixed Effects:
##          (Intr) age    ratio  sexM   bdTmp1 ambTm1
## age      -0.220                                   
## ratio    -0.100 -0.213                            
## sexM     -0.127  0.006 -0.240                     
## bodTemp1 -0.851 -0.266  0.008  0.146              
## ambTemp1  0.139  0.190 -0.082  0.078 -0.356       
## timeOut1 -0.039 -0.100  0.094 -0.031  0.062 -0.206

Plot ratio model

sat1RatioMdl.plot <- plot_model(sat1Ratio.mdl,
           title = "Saturation 1 ratio model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
sat1RatioMdl.plot

Sat1 ~ mass model

sat1Mass.mdl <- lmer(sat1 ~ age + mass + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = mc, REML = FALSE)
summary(sat1Mass.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ age + mass + sex + bodTemp1 + ambTemp1 + timeOut1 + (1 |  
##     family)
##    Data: mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    806.4    829.9   -394.2    788.4       91 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6188 -0.5994  0.1448  0.5711  1.9641 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 102.29   10.114  
##  Residual              97.53    9.876  
## Number of obs: 100, groups:  family, 36
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 96.38436   73.84681   1.305
## age         -3.48163    1.49686  -2.326
## mass         0.05982    0.03622   1.651
## sexM        -4.18234    2.59786  -1.610
## bodTemp1     1.55714    1.88466   0.826
## ambTemp1     1.48949    0.58264   2.556
## timeOut1     0.00317    0.14519   0.022
## 
## Correlation of Fixed Effects:
##          (Intr) age    mass   sexM   bdTmp1 ambTm1
## age      -0.233                                   
## mass     -0.097 -0.164                            
## sexM     -0.111  0.015 -0.355                     
## bodTemp1 -0.849 -0.267  0.008  0.140              
## ambTemp1  0.130  0.176 -0.018  0.062 -0.354       
## timeOut1 -0.042 -0.100  0.120 -0.050  0.063 -0.200

Plot mass model

sat1MassMdl.plot <- plot_model(sat1Mass.mdl,
           title = "Saturation 1 mass model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
sat1MassMdl.plot

Sat1 ~ cube mass model

sat1CubeMass.mdl <- lmer(sat1 ~ age + cubeMass + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = mc, REML = FALSE)
summary(sat1CubeMass.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ age + cubeMass + sex + bodTemp1 + ambTemp1 + timeOut1 +  
##     (1 | family)
##    Data: mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    806.2    829.7   -394.1    788.2       91 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6156 -0.5987  0.1421  0.5741  1.9641 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 100.89   10.044  
##  Residual              97.72    9.886  
## Number of obs: 100, groups:  family, 36
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 52.329008  80.194046   0.653
## age         -3.465349   1.487939  -2.329
## cubeMass     9.462185   5.474152   1.729
## sexM        -4.259462   2.598636  -1.639
## bodTemp1     1.500293   1.884291   0.796
## ambTemp1     1.503016   0.580330   2.590
## timeOut1     0.005269   0.144882   0.036
## 
## Correlation of Fixed Effects:
##          (Intr) age    cubMss sexM   bdTmp1 ambTm1
## age      -0.164                                   
## cubeMass -0.403 -0.154                            
## sexM      0.009  0.012 -0.356                     
## bodTemp1 -0.778 -0.266 -0.008  0.145              
## ambTemp1  0.123  0.176 -0.008  0.059 -0.355       
## timeOut1 -0.078 -0.099  0.125 -0.052  0.060 -0.199

Plot cube mass model

sat1CubeMassMdl.plot <- plot_model(sat1CubeMass.mdl,
           title = "Saturation 1 cube mass model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
sat1CubeMassMdl.plot

Sat1 ~ cube ratio model

sat1CubeRatio.mdl <- lmer(sat1 ~ age + cubeRatio + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = mc, REML = FALSE)
summary(sat1CubeRatio.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ age + cubeRatio + sex + bodTemp1 + ambTemp1 + timeOut1 +  
##     (1 | family)
##    Data: mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    806.2    829.7   -394.1    788.2       91 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.8199 -0.5831  0.1268  0.5632  2.2226 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 120.82   10.992  
##  Residual              91.71    9.576  
## Number of obs: 100, groups:  family, 36
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 55.83330   81.10808   0.688
## age         -3.70334    1.58647  -2.334
## cubeRatio    5.55821    3.57531   1.555
## sexM        -2.00929    2.40637  -0.835
## bodTemp1     1.63945    1.85847   0.882
## ambTemp1     1.27632    0.61714   2.068
## timeOut1    -0.01771    0.14746  -0.120
## 
## Correlation of Fixed Effects:
##           (Intr) age    cubeRt sexM   bdTmp1 ambTm1
## age       -0.177                                   
## cubeRatio -0.413 -0.210                            
## sexM      -0.202 -0.074  0.152                     
## bodTemp1  -0.763 -0.248  0.002  0.152              
## ambTemp1   0.163  0.194 -0.180  0.024 -0.329       
## timeOut1  -0.036 -0.077  0.003  0.002  0.070 -0.186

Plot cube mass model

sat1CubeRatioMdl.plot <- plot_model(sat1CubeRatio.mdl,
           title = "Saturation 1 cube ratio model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
sat1CubeRatioMdl.plot

### Now let’s compare the fit of models with different condition indices

indices.mdls <- c(sat1Residual.mdl, sat1Ratio.mdl, sat1Mass.mdl, sat1CubeMass.mdl, sat1CubeRatio.mdl)
indices.modnames <- c("Residual", "Ratio", "Mass", "Cube Mass", "Cube Ratio")
confset(indices.mdls, modnames = indices.modnames, method = "ordinal")

## 
## Confidence set for the best model
## 
## Method:   ordinal ranking based on delta AIC
## 
## Models with substantial weight:
##            K   AICc Delta_AICc AICcWt
## Ratio      9 807.32       0.00   0.26
## Residual   9 807.43       0.11   0.25
## Cube Mass  9 808.22       0.89   0.17
## Cube Ratio 9 808.22       0.90   0.17
## Mass       9 808.42       1.09   0.15
## 
## 
## Models with some weight:
##      K AICc Delta_AICc AICcWt
## 
## 
## Models with little weight:
##      K AICc Delta_AICc AICcWt
## 
## 
## Models with no weight:
##      K AICc Delta_AICc AICcWt

Using AICc Weight (AICcWt), the confset() function calculates the probability of each model given the data and the other candidate models. So in this case, we have more confidence in the ratio model’s ability to explain the variation in the data, even though AICc values are only slightly different.

HUE 1 MODELS

hue1.mdl <- lmer(hue1*100 ~ age + ratio + sex + bodTemp1 + ambTemp1 + (1|family), 
                    data = mc, REML = FALSE)

hue1Mdl.plot <- plot_model(hue1.mdl,
           title = "hue 1 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
hue1Mdl.plot

TIME 2 DATA FRAME

#create a time 2 data frame (N=103)
time2.mc <- mc %>% select(row:sex, ends_with("2"), ambFinal, broodSize, date:skull, resids, ratio)

#I'm dropping ambTemp2 from data frame before dropping NAs. 
#ambTemp2 has many more NAs. We can use ambFinal instead.
#This boosts sample size from 32 to 99 obs. 
time2.mc <- time2.mc %>% select(!ambTemp2) 

#Now drop NAs after removing ambTemp2
time2.mc <- drop_na(time2.mc)

At time 1, we have 78 females and 55 males. At time 2, we have 64 females and 39 males.

SATURATION 2 MODELS

sat2.mdl <- lmer(sat2 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 + (1|family), 
                    data = time2.mc, REML = FALSE)
summary(sat2.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat2 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 + (1 |  
##     family)
##    Data: time2.mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    834.5    858.2   -408.2    816.5       94 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.81870 -0.61972  0.01104  0.52422  2.44593 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 144.83   12.034  
##  Residual              89.75    9.474  
## Number of obs: 103, groups:  family, 38
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 84.243482  66.498688   1.267
## age         -3.949011   1.631450  -2.421
## ratio        1.828675   2.496676   0.732
## sexM         4.049973   2.466390   1.642
## bodTemp2     2.230622   1.679602   1.328
## ambFinal     1.806990   0.567272   3.185
## timeOut2     0.005522   0.064403   0.086
## 
## Correlation of Fixed Effects:
##          (Intr) age    ratio  sexM   bdTmp2 ambFnl
## age      -0.336                                   
## ratio    -0.004 -0.124                            
## sexM     -0.021  0.038 -0.378                     
## bodTemp2 -0.778 -0.258 -0.168  0.064              
## ambFinal  0.047  0.049  0.095  0.083 -0.257       
## timeOut2 -0.101 -0.037  0.137 -0.129  0.057 -0.097

Plot sat2 model

sat2Mdl.plot <- plot_model(sat2.mdl,
           title = "Saturation 2 ratio model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
sat2Mdl.plot

HUE 2 MODELS

hue2.mdl <- lmer(hue2*100 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 + (1|family), 
                    data = time2.mc, REML = FALSE)
summary(hue2.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: hue2 * 100 ~ age + ratio + sex + bodTemp2 + ambFinal + timeOut2 +  
##     (1 | family)
##    Data: time2.mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    565.2    588.9   -273.6    547.2       94 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.47645 -0.47349  0.00881  0.53622  2.39858 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 6.795    2.607   
##  Residual             7.731    2.780   
## Number of obs: 103, groups:  family, 38
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept) 114.52826   17.46943   6.556
## age           0.00711    0.38966   0.018
## ratio        -1.63567    0.66936  -2.444
## sexM          0.03974    0.69593   0.057
## bodTemp2     -0.87323    0.46688  -1.870
## ambFinal     -0.27935    0.13579  -2.057
## timeOut2     -0.04941    0.01737  -2.845
## 
## Correlation of Fixed Effects:
##          (Intr) age    ratio  sexM   bdTmp2 ambFnl
## age      -0.210                                   
## ratio    -0.001 -0.144                            
## sexM     -0.022  0.033 -0.338                     
## bodTemp2 -0.826 -0.306 -0.161  0.055              
## ambFinal  0.097  0.080  0.101  0.107 -0.290       
## timeOut2 -0.115 -0.056  0.094 -0.115  0.084 -0.124

Plot hue 2 model

hue2Mdl.plot <- plot_model(hue2.mdl,
           title = "hue 2 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
hue2Mdl.plot

Wrap sat and hue model plots

ratioMdlWrap <- ggarrange(sat1RatioMdl.plot, hue1Mdl.plot, sat2Mdl.plot, hue2Mdl.plot)
ratioMdlWrap

#ggsave("ratioMdlWrap.png", plot = ratioMdlWrap, device = "png", dpi = 300, width = 10, height = 8)

Birds with both time1 and time2 data ONLY

The semi_join() function returns all rows from mc where there are matching IDs in time2.mc, keeping just columns from mc.

bothTimes.mc <- semi_join(mc,time2.mc, by="id") #N=103

BOTH TIMES SAT1 RATIO MODEL

sat1.BT.mdl <- lmer(sat1 ~ age + ratio + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = bothTimes.mc, REML = FALSE)
summary(sat1.BT.mdl)

## Linear mixed model fit by maximum likelihood  ['lmerMod']
## Formula: sat1 ~ age + ratio + sex + bodTemp1 + ambTemp1 + timeOut1 + (1 |  
##     family)
##    Data: bothTimes.mc
## 
##      AIC      BIC   logLik deviance df.resid 
##    715.8    738.2   -348.9    697.8       80 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.0229 -0.5839  0.1569  0.6019  1.8685 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  family   (Intercept) 122.56   11.071  
##  Residual              85.07    9.223  
## Number of obs: 89, groups:  family, 33
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept) 86.87942   77.80080   1.117
## age         -3.81232    1.60012  -2.383
## ratio        3.86820    2.67839   1.444
## sexM        -3.62263    2.51583  -1.440
## bodTemp1     2.05087    1.95853   1.047
## ambTemp1     1.26256    0.63269   1.996
## timeOut1    -0.01961    0.15890  -0.123
## 
## Correlation of Fixed Effects:
##          (Intr) age    ratio  sexM   bdTmp1 ambTm1
## age      -0.245                                   
## ratio    -0.104 -0.231                            
## sexM     -0.143  0.032 -0.259                     
## bodTemp1 -0.847 -0.249  0.028  0.152              
## ambTemp1  0.106  0.187 -0.101  0.099 -0.330       
## timeOut1  0.006 -0.099  0.053 -0.055  0.030 -0.203

Plot both times sat1 model

sat1.BT.mdl.plot <- plot_model(sat1.BT.mdl,
           title = "Both times sat1 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
sat1.BT.mdl.plot

Both times hue1 ratio model

hue1.BT.mdl <- lmer(hue1*100 ~ age + ratio + sex + bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = bothTimes.mc, REML = FALSE)

Plot both times hue1 model

hue1.BT.mdl.plot <- plot_model(hue1.BT.mdl,
           title = "Both times hue1 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)
hue1.BT.mdl.plot

Wrap both times model plots

bothTimes.mdlWrap <- ggarrange(sat1.BT.mdl.plot, hue1.BT.mdl.plot, sat2Mdl.plot, hue2Mdl.plot)
bothTimes.mdlWrap

#ggsave("bothTimesMdlWrap.png", plot = bothTimes.mdlWrap, device = "png", dpi = 300, width = 10, height = 8)

WHAT RESPONSE VARIABLES TO USE?

Plot saturation by hue

mc %>% ggplot(aes(x=sat1, y=hue1))+
  geom_jitter()+
  geom_smooth()

## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

Separate yellow and magenta models for time 1

Mean yellow 1 model

meanY.mdl <- lmer(meanY1 ~ age + ratio + sex + 
                        bodTemp1 + ambTemp1 + (1|family), 
                    data = mc, REML = FALSE)

Plot mean yellow 1 model

meanY.plot <- plot_model(meanY.mdl,
           title = "mean yellow 1 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Mean magenta 1 model

meanM.mdl <- lmer(meanM1 ~ age + ratio + sex + 
                        bodTemp1 + ambTemp1 + (1|family), 
                    data = mc, REML = FALSE)

Plot mean magenta 1 model

meanM.plot <- plot_model(meanM.mdl,
           title = "mean magenta 1 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Wrap plot of mean Y1 & mean M1 models

ggarrange(meanY.plot,meanM.plot)

Separate yellow and magenta models for time 2

Mean yellow 2 model

meanY2.mdl <- lmer(meanY2 ~ age + ratio + sex + 
                        bodTemp2 + ambFinal + timeOut2 + (1|family), 
                    data = time2.mc, REML = FALSE)

Plot mean yellow 2 model

meanY2.plot <- plot_model(meanY2.mdl,
           title = "mean yellow 2 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Mean magenta 2 model

meanM2.mdl <- lmer(meanM2 ~ age + ratio + sex + 
                        bodTemp2 + ambFinal + timeOut2 + (1|family), 
                    data = time2.mc, REML = FALSE)

Plot mean magenta 2 model

meanM2.plot <- plot_model(meanM2.mdl,
           title = "mean magenta 2 model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Wrap plots of mean Y2 & M2 models

ggarrange(meanY2.plot,meanM2.plot)

Wrap model plots for mean colors and hue at time1 and time2

meanColorMdlWrap <- ggarrange(meanY.plot,meanM.plot,hue1Mdl.plot, meanY2.plot,meanM2.plot,hue2Mdl.plot)
meanColorMdlWrap

#ggsave("meanColorMdlWrap.png", plot = meanColorMdlWrap, device = "png", dpi = 300, width = 10, height = 8)

Include older birds at time 1

Data frame including birds 23-33 days old

older.mc <- allAges.mc %>% filter(between(age,23,33))
#create mass/tarsus ratio column
older.mc$ratio <- older.mc$mass/older.mc$tarsus

Mean yellow 1 model with older birds

olderMeanY.mdl <- lmer(meanY1 ~ age + ratio + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = older.mc, REML = FALSE)

Plot mean yellow 1 model with older birds

olderMeanY.plot <- plot_model(olderMeanY.mdl,
           title = "mean yellow model with older birds",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Mean magenta 1 model with older birds

olderMeanM.mdl <- lmer(meanM1 ~ age + ratio + sex + 
                        bodTemp1 + ambTemp1 + (1|family), 
                    data = older.mc, REML = FALSE)

Plot mean magenta 1 model with older birds

olderMeanM.plot <- plot_model(olderMeanM.mdl,
           title = "mean magenta model with older birds",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

Wrap mean colors with older birds models

ggarrange(olderMeanY.plot,olderMeanM.plot,meanY.plot,meanM.plot)

Create age bins factor for older bird data frame

#define bin breaks
binBreaks <- c(22,25,28,33) #define breaks
#name bins
bins <- c("bin1","bin2","bin3") #label bins
#assign age variable to bins
older.mc$ageBins <- cut(older.mc$age,breaks = binBreaks,labels = bins) #vectorize
#create a new column for age bins
older.mc$ageBins <- as.factor(older.mc$ageBins)
#checkBins <- data.frame(older.mc$id, older.mc$age, older.mc$ageBins)

ageBinsMeanY.mdl <- lmer(meanY1 ~ ageBins + ratio + sex + 
                        bodTemp1 + ambTemp1 + timeOut1 + (1|family), 
                    data = older.mc, REML = FALSE)

ageBinsMeanY.plot <- plot_model(ageBinsMeanY.mdl,
           title = "age bin mean yellow model",
           bpe = "median",
           vline.color = "black",
           show.values = TRUE, 
           value.offset = .3,
           line.size = 1.2)

ggarrange(olderMeanY.plot, ageBinsMeanY.plot)

Birds in the oldest category (29-33) are driving the negative age effect.

Comparison of model fit for continuous age and age bin models

ageBins.mdls <- c(olderMeanY.mdl, ageBinsMeanY.mdl)
ageBins.modnames <- c("continuous age", "age bins")
confset(ageBins.mdls, modnames = ageBins.modnames, method = "ordinal")

## 
## Confidence set for the best model
## 
## Method:   ordinal ranking based on delta AIC
## 
## Models with substantial weight:
##                K    AICc Delta_AICc AICcWt
## continuous age 9 1021.89          0   0.95
## 
## 
## Models with some weight:
##           K    AICc Delta_AICc AICcWt
## age bins 10 1027.84       5.95   0.05
## 
## 
## Models with little weight:
##      K AICc Delta_AICc AICcWt
## 
## 
## Models with no weight:
##      K AICc Delta_AICc AICcWt

There’s no difference in the model probabilities.